Diversity Analysis
Diversity Analysis
Import libraries and data
library(phyloseq)
library(ape)
library(tidyverse)## Warning: package 'ggplot2' was built under R version 4.3.2## Warning: package 'tidyr' was built under R version 4.3.2## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::where() masks ape::where()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(picante)## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.3.2## This is vegan 2.6-4
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapselibrary(ARTool)
library(ggsignif)
library(pairwiseAdonis)## Loading required package: clusterlibrary(ggpubr)##
## Attaching package: 'ggpubr'
##
## The following object is masked from 'package:ape':
##
## rotatestarch <- read_rds("../21_phyloseq/starch_rarefied_phyloseq.RDS")Separate human and mouse for the independent analysis
human <- subset_samples(starch, cohort == "human")
mouse <- subset_samples(starch, cohort == "mouse")pal <- c("#DD4124FF", "#faa478", "#6ac1de", "#0269a1FF")
pal2 <- c("#00346EFF", "#007CBFFF", "#06ABDFFF", "#EDD03EFF", "#F5A200FF", "#6D8600FF", "#424D0CFF")
pal3 <- c("#FDAE61FF", "#ABDDA4FF")Alpha diversity
richness <- estimate_richness(starch)
rownames(richness) <- rownames(starch@sam_data)starch_rich <- merge(starch@sam_data, richness, by = 0) %>% dplyr::rename(sample_id = Row.names)WIP:
otu <- otu_table(starch) %>% as.matrix() data_evenness <- diversity(otu) / log(specnumber(otu))Shannon and Chao1
Statistical Analysis
Merged
shannon_kw <- kruskal.test(Shannon ~ group, data = starch_rich)
shannon_kw##
## Kruskal-Wallis rank sum test
##
## data: Shannon by group
## Kruskal-Wallis chi-squared = 26.916, df = 3, p-value = 6.132e-06if (shannon_kw$p.value < 0.05){
pairwise.wilcox.test(starch_rich$Shannon, starch_rich$group,
p.adjust = "bonferroni")
}##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: starch_rich$Shannon and starch_rich$group
##
## human-C human-RS mouse-C
## human-RS 0.058 - -
## mouse-C 6.1e-06 0.019 -
## mouse-RS 9.8e-05 0.218 1.000
##
## P value adjustment method: bonferroniThere is a significant difference (p-val < 0.05) in the Shannon index between: - mouse-C & human-C - mouse-C & human-RS - mouse-RS & human-C
shannon_kw <- kruskal.test(Shannon ~ gender, data = starch_rich)
shannon_kw##
## Kruskal-Wallis rank sum test
##
## data: Shannon by gender
## Kruskal-Wallis chi-squared = 5.6784, df = 1, p-value = 0.01718Overall shows a significant difference in the Shannon index between genders.
chao1_kw <- kruskal.test(Chao1 ~ group, data = starch_rich)
chao1_kw##
## Kruskal-Wallis rank sum test
##
## data: Chao1 by group
## Kruskal-Wallis chi-squared = 10.757, df = 3, p-value = 0.01311if (chao1_kw$p.value < 0.05){
pairwise.wilcox.test(starch_rich$Chao1, starch_rich$group,
p.adjust = "bonferroni")
}## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: starch_rich$Chao1 and starch_rich$group
##
## human-C human-RS mouse-C
## human-RS 1.000 - -
## mouse-C 0.018 0.770 -
## mouse-RS 0.045 1.000 1.000
##
## P value adjustment method: bonferroniDifference in Chao1 index between: human-c & mouse-c and human-C & mouse-rs No difference between: human-RS and mouse-RS!!
chao1_kw <- kruskal.test(Chao1 ~ gender, data = starch_rich)
chao1_kw##
## Kruskal-Wallis rank sum test
##
## data: Chao1 by gender
## Kruskal-Wallis chi-squared = 3.33, df = 1, p-value = 0.06803if (chao1_kw$p.value < 0.05){
pairwise.wilcox.test(starch_rich$Chao1, starch_rich$gender,
p.adjust = "bonferroni")
}No significant difference between genders.
Mouse
starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(Shannon ~ trt, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Shannon by trt
## Kruskal-Wallis chi-squared = 1.2027, df = 1, p-value = 0.2728No significant difference in Shannon index between treatments.
starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(Chao1 ~ trt, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Chao1 by trt
## Kruskal-Wallis chi-squared = 0.98674, df = 1, p-value = 0.3205No significant difference in Chao1 index between treatments.
starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(Shannon ~ gender, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Shannon by gender
## Kruskal-Wallis chi-squared = 8.3409, df = 1, p-value = 0.003876Significant difference in Shannon index between genders.
starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(Chao1 ~ gender, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Chao1 by gender
## Kruskal-Wallis chi-squared = 5.6614, df = 1, p-value = 0.01734Significant difference in Chao1 index between genders.
starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(Shannon ~ starch_type, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Shannon by starch_type
## Kruskal-Wallis chi-squared = 8.2792, df = 4, p-value = 0.08187No significant difference in Shannon index between starch types
mouse_kw <- starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(Chao1 ~ starch_type, data = .)mouse_rich <- starch_rich %>%
filter(cohort == "mouse")
chao1_kw <- mouse_rich %>%
kruskal.test(Chao1 ~ starch_type, data = .)
if (chao1_kw$p.value < 0.05){
pairwise.wilcox.test(mouse_rich$Chao1, mouse_rich$starch_type,
p.adjust = "bonferroni")
}## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties
## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: mouse_rich$Chao1 and mouse_rich$starch_type
##
## BEP CKP CTL LEN
## CKP 0.068 - - -
## CTL 1.000 0.092 - -
## LEN 1.000 0.459 1.000 -
## PTB 1.000 0.028 1.000 1.000
##
## P value adjustment method: bonferroniSignificant difference in mean Chao1 index between CKP and PTB.
Human
starch_rich %>%
filter(cohort == "human") %>%
kruskal.test(Shannon ~ trt, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Shannon by trt
## Kruskal-Wallis chi-squared = 6.5108, df = 1, p-value = 0.01072Significant difference in Shannon index between treatments
starch_rich %>%
filter(cohort == "human") %>%
kruskal.test(Chao1 ~ trt, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Chao1 by trt
## Kruskal-Wallis chi-squared = 1.5583, df = 1, p-value = 0.2119No significant difference in Chao1 index between treatments.
starch_rich %>%
filter(cohort == "human") %>%
kruskal.test(Shannon ~ gender, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Shannon by gender
## Kruskal-Wallis chi-squared = 4.6282, df = 1, p-value = 0.03145Significant difference in Shannon index between genders.
starch_rich %>%
filter(cohort == "human") %>%
kruskal.test(Chao1 ~ gender, data = .)##
## Kruskal-Wallis rank sum test
##
## data: Chao1 by gender
## Kruskal-Wallis chi-squared = 1.3475, df = 1, p-value = 0.2457No significant difference in Chao1 index between genders.
Plots
my_comparisons <- list(c("human-C", "human-RS"),
c("mouse-C", "mouse-RS"),
c("mouse-RS","human-RS"),
c("mouse-C","human-C"))
starch_comparisons <- list(c("CTL", "CKP"),
c("CKP", "LEN"),
c("CKP", "BEP"),
c("CKP", "PTB")
)Merged
w <- 2400
h <- 2000starch_rich %>%
ggboxplot(x = "group", y = "Shannon", color = "group",
palette = pal) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank()
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_group.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
ggboxplot(x = "gender", y = "Shannon", color = "gender",
palette = pal3) +
stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank()
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
ggboxplot(x = "gender", y = "Shannon", color = "gender",
palette = pal3, facet.by = "group") +
stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank()
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender_facet.png", units = "px", width = w, height = h, create.dir = TRUE)starch_rich %>%
ggboxplot(x = "group", y = "Chao1", color = "group",
palette = pal) +
stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank()
)## Warning in wilcox.test.default(c(94.3333333333333, 130.5, 83, 124, 85, 97, :
## cannot compute exact p-value with ties## Warning in wilcox.test.default(c(67, 75, 58, 71, 69.4285714285714, 72, 69.5, :
## cannot compute exact p-value with ties
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_group.png", units = "px", width = w, height = h, create.dir = TRUE)## Warning in wilcox.test.default(c(94.3333333333333, 130.5, 83, 124, 85, 97, :
## cannot compute exact p-value with ties
## Warning in wilcox.test.default(c(94.3333333333333, 130.5, 83, 124, 85, 97, :
## cannot compute exact p-value with tiesstarch_rich %>%
ggboxplot(x = "gender", y = "Chao1", color = "gender",
palette = pal3) +
stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank()
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
ggboxplot(x = "gender", y = "Chao1", color = "gender",
palette = pal3, facet.by = "group") +
stat_compare_means(label = "p.signif", label.y = 4.2, label.x = 1.5) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank()
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender_facet.png", units = "px", width = w, height = h, create.dir = TRUE)Mouse
starch_rich %>%
filter(cohort == "mouse") %>%
ggboxplot(x = "gender", y = "Shannon", color = "gender",
palette = pal3) +
facet_wrap(.~trt) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
filter(cohort == "mouse") %>%
ggboxplot(x = "trt", y = "Shannon", color = "trt",
palette = rev(pal[3:4])) +
facet_wrap(.~gender) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
filter(cohort == "mouse") %>%
ggboxplot(x = "starch_type", y = "Shannon", color = "starch_type",
palette = pal2) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
xlab("Starch Type") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)starch_rich %>%
filter(cohort == "mouse") %>%
ggboxplot(x = "gender", y = "Chao1", color = "gender",
palette = pal3) +
facet_wrap(.~trt) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
filter(cohort == "mouse") %>%
ggboxplot(x = "trt", y = "Chao1", color = "trt",
palette = rev(pal[3:4])) +
facet_wrap(.~gender) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
filter(cohort == "mouse") %>%
ggboxplot(x = "starch_type", y = "Chao1", color = "starch_type",
palette = pal2) +
stat_compare_means(comparisons = starch_comparisons,
method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
xlab("Starch Type") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none"
)## Warning: Computation failed in `stat_signif()`.
## Caused by error in `kruskal.test.default()`:
## ! 'x' and 'g' must have the same length
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)## Warning: Computation failed in `stat_signif()`.
## Caused by error in `kruskal.test.default()`:
## ! 'x' and 'g' must have the same lengthHuman
starch_rich %>%
filter(cohort == "human") %>%
ggboxplot(x = "gender", y = "Shannon", color = "gender",
palette = pal3) +
facet_wrap(.~trt) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
filter(cohort == "human") %>%
ggboxplot(x = "trt", y = "Shannon", color = "trt",
palette = pal) +
facet_wrap(.~gender) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/shannon_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)starch_rich %>%
filter(cohort == "human") %>%
ggboxplot(x = "gender", y = "Chao1", color = "gender",
palette = pal3) +
facet_wrap(.~trt) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
filter(cohort == "human") %>%
ggboxplot(x = "trt", y = "Chao1", color = "trt",
palette = pal) +
facet_wrap(.~gender) +
stat_compare_means(method = "kruskal.test",
hide.ns = TRUE,
label = "p.signif") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "right"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/chao1_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)Faith’s Phylogenetic Diversity
phylo_dist <- pd(t(otu_table(starch)), phy_tree(starch),
include.root=F)
# add PD to metadata table
starch_rich$PD <- phylo_dist$PDStatistical Analysis
Group
pd_kw <- kruskal.test(PD ~ group, data = starch_rich)
pd_kw##
## Kruskal-Wallis rank sum test
##
## data: PD by group
## Kruskal-Wallis chi-squared = 4.5904, df = 3, p-value = 0.2044if (pd_kw$p.value < 0.05){
pairwise.wilcox.test(starch_rich$PD, starch_rich$group,
p.adjust = "bonferroni")
}There is no significant difference in Faith’s phylogenetic distance between the groups.
Gender
kruskal.test(PD ~ gender, data = starch_rich)##
## Kruskal-Wallis rank sum test
##
## data: PD by gender
## Kruskal-Wallis chi-squared = 1.2316, df = 1, p-value = 0.2671No difference between genders.
Mouse
starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(PD ~ gender, data = .)##
## Kruskal-Wallis rank sum test
##
## data: PD by gender
## Kruskal-Wallis chi-squared = 3.1929, df = 1, p-value = 0.07396No difference between genders.
starch_rich %>%
filter(cohort == "mouse") %>%
kruskal.test(PD ~ starch_type, data = .)##
## Kruskal-Wallis rank sum test
##
## data: PD by starch_type
## Kruskal-Wallis chi-squared = 6.8053, df = 4, p-value = 0.1465No difference between starch types.
Human
starch_rich %>%
filter(cohort == "human") %>%
kruskal.test(PD ~ gender, data = .)##
## Kruskal-Wallis rank sum test
##
## data: PD by gender
## Kruskal-Wallis chi-squared = 1.0385, df = 1, p-value = 0.3082No difference between genders.
Plots
Merged
starch_rich %>%
ggboxplot(x = "group", y = "PD", color = "group",
palette = pal) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none",
axis.title.x = element_blank()
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/pd_group.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
ggboxplot(x = "gender", y = "PD", color = "gender",
palette = pal3) +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/pd_gender.png", units = "px", width = w, height = h, create.dir = TRUE)
starch_rich %>%
ggboxplot(x = "gender", y = "PD", color = "gender",
palette = pal3, facet.by = "group") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = "none"
)
ggsave("../30_diversity-analysis/plots/alpha-diversity/pd_gender_facet.png", units = "px", width = w, height = h, create.dir = TRUE)Beta Diversity
starchdf <- as.data.frame(sample_data(starch))
humandf <- as.data.frame(sample_data(human))
mousedf <- as.data.frame(sample_data(mouse))Bray Curtis
Merged
bc_dm <- phyloseq::distance(starch, method="bray")
adonis2(bc_dm ~ group*gender, data = data.frame(starchdf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bc_dm ~ group * gender, data = data.frame(starchdf))
## Df SumOfSqs R2 F Pr(>F)
## group 3 9.325 0.25341 13.304 0.001 ***
## gender 1 3.421 0.09296 14.641 0.001 ***
## group:gender 3 2.323 0.06313 3.314 0.001 ***
## Residual 93 21.729 0.59050
## Total 100 36.797 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(bc_dm ~ trt*cohort*gender, data = data.frame(starchdf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bc_dm ~ trt * cohort * gender, data = data.frame(starchdf))
## Df SumOfSqs R2 F Pr(>F)
## trt 1 1.737 0.04722 7.4361 0.001 ***
## cohort 1 7.253 0.19710 31.0412 0.001 ***
## gender 1 3.451 0.09379 14.7715 0.001 ***
## trt:cohort 1 0.304 0.00827 1.3028 0.208
## trt:gender 1 0.983 0.02670 4.2055 0.001 ***
## cohort:gender 1 0.956 0.02597 4.0898 0.001 ***
## trt:cohort:gender 1 0.385 0.01046 1.6468 0.109
## Residual 93 21.729 0.59050
## Total 100 36.797 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant difference between treatments, cohorts, genders, at least two trt:gender pair, and at least two cohort:gender pair.
pcoa_bc <- ordinate(starch, method="PCoA", distance=bc_dm)
plot_ordination(starch, pcoa_bc, color = "group") +
scale_color_manual(values = pal) +
labs(col = "Group") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_group.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(starch, pcoa_bc, color = "gender") +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender.png", units = "px", width = w, height = h, create.dir = TRUE)Human
human_bc_dm <- distance(human, method="bray")
adonis2(human_bc_dm ~ trt*gender, data = data.frame(humandf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = human_bc_dm ~ trt * gender, data = data.frame(humandf))
## Df SumOfSqs R2 F Pr(>F)
## trt 1 0.2863 0.04489 1.0845 0.322
## gender 1 0.4024 0.06310 1.5246 0.028 *
## trt:gender 1 0.1455 0.02281 0.5512 0.993
## Residual 21 5.5430 0.86919
## Total 24 6.3772 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is only a significant difference in community structure between
genders. I should probably consider here the effect of
trtseq (fixed effect) or individuals as random error….
human_pcoa_bc <- ordinate(human, method="PCoA", distance=human_bc_dm)
plot_ordination(human, human_pcoa_bc, color = "trt") +
scale_color_manual(values = pal) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(human, human_pcoa_bc, color = "gender") +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(human, human_pcoa_bc, color = "gender") +
scale_color_manual(values = pal3) +
facet_grid(.~trt) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender_facet_human.png", units = "px", width = w, height = h, create.dir = TRUE)Mouse
mouse_bc_dm <- distance(mouse, method="bray")
adonis2(mouse_bc_dm ~ trt*starch_type*gender, data = data.frame(mousedf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = mouse_bc_dm ~ trt * starch_type * gender, data = data.frame(mousedf))
## Df SumOfSqs R2 F Pr(>F)
## trt 1 0.4224 0.01937 2.6609 0.030 *
## starch_type 3 3.0381 0.13934 6.3789 0.001 ***
## gender 1 4.2557 0.19518 26.8060 0.001 ***
## trt:gender 1 0.8667 0.03975 5.4595 0.001 ***
## starch_type:gender 3 2.7427 0.12579 5.7587 0.001 ***
## Residual 66 10.4781 0.48056
## Total 75 21.8038 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pw <- pairwiseAdonis::pairwise.adonis2(mouse_bc_dm ~ trt*starch_type*gender,
data = data.frame(mousedf))
print("significantly different pairs: ")## [1] "significantly different pairs: "for (i in range(2,(length(pw)))){
pair <- pw[[i]]
if (pair$`Pr(>F)`[1] < 0.05) {
print(names(pw)[i])
}
}## [1] "RS_vs_C"
## [1] "RS_vs_C"There is a significant difference in community structure between treatments, at least 2 starch types, genders, at least two trt:gender pairs, and at least two starch_type:gender pairs.
mouse_pcoa_bc <- ordinate(mouse, method="PCoA", distance=mouse_bc_dm)
plot_ordination(mouse, mouse_pcoa_bc, color = "trt") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal[3:4]) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_bc, color = "starch_type") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal2) +
labs(col = "Starch Type") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_bc, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Starch Type") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/bc_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)Jaccard
Merged
jc_dm <- phyloseq::distance(starch, method = "jaccard", binary = TRUE)
adonis2(jc_dm ~ group*gender, data = data.frame(starchdf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = jc_dm ~ group * gender, data = data.frame(starchdf))
## Df SumOfSqs R2 F Pr(>F)
## group 3 7.162 0.17428 7.0604 0.001 ***
## gender 1 1.194 0.02906 3.5315 0.001 ***
## group:gender 3 1.292 0.03144 1.2736 0.042 *
## Residual 93 31.444 0.76522
## Total 100 41.092 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(jc_dm ~ cohort*trt*gender, data = data.frame(starchdf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = jc_dm ~ cohort * trt * gender, data = data.frame(starchdf))
## Df SumOfSqs R2 F Pr(>F)
## cohort 1 6.352 0.15458 18.7866 0.001 ***
## trt 1 0.439 0.01070 1.2999 0.094 .
## gender 1 1.204 0.02930 3.5609 0.001 ***
## cohort:trt 1 0.360 0.00877 1.0655 0.260
## cohort:gender 1 0.641 0.01559 1.8947 0.010 **
## trt:gender 1 0.375 0.00912 1.1082 0.237
## cohort:trt:gender 1 0.277 0.00673 0.8180 0.839
## Residual 93 31.444 0.76522
## Total 100 41.092 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant difference in community composition between cohorts, genders, and at least 2 cohort:gender pairs.
pcoa_jc <- ordinate(starch, method="PCoA", distance=jc_dm)
plot_ordination(starch, pcoa_jc, color = "group") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal) +
labs(col = "Group") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/jc_group.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(starch, pcoa_jc, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/jc_gender.png", units = "px", width = w, height = h, create.dir = TRUE)Human
human_jc_dm <- distance(human, method = "jaccard", binary = TRUE)
adonis2(human_jc_dm ~ trt*gender, data = data.frame(humandf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = human_jc_dm ~ trt * gender, data = data.frame(humandf))
## Df SumOfSqs R2 F Pr(>F)
## trt 1 0.2831 0.03702 0.8769 0.762
## gender 1 0.4030 0.05271 1.2484 0.092 .
## trt:gender 1 0.1805 0.02361 0.5592 1.000
## Residual 21 6.7794 0.88666
## Total 24 7.6461 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Community composition is not significantly different between treatments or gender.
human_pcoa_jc <- ordinate(human, method="PCoA", distance=human_jc_dm)
plot_ordination(human, human_pcoa_jc, color = "trt") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/jc_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(human, human_pcoa_jc, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/jc_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)Mouse
mouse_jc_dm <- distance(mouse, method = "jaccard", binary = TRUE)
adonis2(mouse_jc_dm ~ gender*trt*starch_type, data = data.frame(mousedf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = mouse_jc_dm ~ gender * trt * starch_type, data = data.frame(mousedf))
## Df SumOfSqs R2 F Pr(>F)
## gender 1 1.4456 0.05335 4.4497 0.001 ***
## trt 1 0.5127 0.01892 1.5781 0.004 **
## starch_type 3 1.7031 0.06286 1.7475 0.001 ***
## gender:trt 1 0.4717 0.01741 1.4518 0.005 **
## gender:starch_type 3 1.5197 0.05609 1.5593 0.001 ***
## Residual 66 21.4413 0.79137
## Total 75 27.0939 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pw <- pairwiseAdonis::pairwise.adonis2(mouse_jc_dm ~ starch_type, data = data.frame(mousedf))
print("significantly different pairs: ")## [1] "significantly different pairs: "for (i in range(2,(length(pw)))){
pair <- pw[[i]]
if (pair$`Pr(>F)`[1] < 0.05) {
print(names(pw)[i])
}
}## [1] "PTB_vs_BEP"
## [1] "CKP_vs_CTL"Community composition is significantly different between genders, treatments, starch types, at least 2 gender:trt pairs and at least 2 gender:starch_type pairs. The significantly different starch_type groups are PTB-BEP and CKP-CTL.
mouse_pcoa_jc <- ordinate(mouse, method="PCoA", distance=mouse_jc_dm)
plot_ordination(mouse, mouse_pcoa_jc, color = "trt") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal[3:4]) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/jc_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_jc, color = "starch_type") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal2) +
labs(col = "Starch Type") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/jc_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_jc, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/jc_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)Unweighted Unifrac
Merged
uu_dm <- phyloseq::distance(starch, method = "unifrac")## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]adonis2(uu_dm ~ group*gender, data = data.frame(starchdf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = uu_dm ~ group * gender, data = data.frame(starchdf))
## Df SumOfSqs R2 F Pr(>F)
## group 3 4.1175 0.21931 9.5810 0.001 ***
## gender 1 0.6622 0.03527 4.6226 0.001 ***
## group:gender 3 0.6730 0.03585 1.5660 0.031 *
## Residual 93 13.3224 0.70958
## Total 100 18.7751 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(uu_dm ~ cohort*trt*gender, data = data.frame(starchdf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = uu_dm ~ cohort * trt * gender, data = data.frame(starchdf))
## Df SumOfSqs R2 F Pr(>F)
## cohort 1 3.8031 0.20256 26.5482 0.001 ***
## trt 1 0.1644 0.00876 1.1478 0.282
## gender 1 0.6630 0.03532 4.6286 0.001 ***
## cohort:trt 1 0.1491 0.00794 1.0409 0.356
## cohort:gender 1 0.3743 0.01994 2.6129 0.008 **
## trt:gender 1 0.1801 0.00959 1.2573 0.225
## cohort:trt:gender 1 0.1186 0.00632 0.8278 0.642
## Residual 93 13.3224 0.70958
## Total 100 18.7751 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant difference in community composition and phylogeny between cohorts, genders and at least 2 group:gender pairs.
pcoa_uu <- ordinate(starch, method="PCoA", distance=uu_dm)
plot_ordination(starch, pcoa_uu, color = "group") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal) +
labs(col = "Group") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/uu_group.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(starch, pcoa_uu, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/uu_gender.png", units = "px", width = w, height = h, create.dir = TRUE)Human
human_uu_dm <- distance(human, method = "unifrac")## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]adonis2(human_uu_dm ~ trt*gender, data = data.frame(humandf), permutations = 9999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = human_uu_dm ~ trt * gender, data = data.frame(humandf), permutations = 9999)
## Df SumOfSqs R2 F Pr(>F)
## trt 1 0.1656 0.04265 1.0377 0.3880
## gender 1 0.2459 0.06336 1.5414 0.0456 *
## trt:gender 1 0.1196 0.03081 0.7497 0.8319
## Residual 21 3.3503 0.86318
## Total 24 3.8814 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant difference in community composition and phylogeny between genders.
human_pcoa_uu <- ordinate(human, method="PCoA", distance=human_uu_dm)
plot_ordination(human, human_pcoa_uu, color = "trt") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/uu_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(human, human_pcoa_uu, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/uu_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)Mouse
mouse_uu_dm <- distance(mouse, method = "unifrac")## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]adonis2(mouse_uu_dm ~ trt*gender*starch_type, data = data.frame(mousedf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = mouse_uu_dm ~ trt * gender * starch_type, data = data.frame(mousedf))
## Df SumOfSqs R2 F Pr(>F)
## trt 1 0.1488 0.01342 1.1175 0.301
## gender 1 0.7906 0.07128 5.9357 0.001 ***
## starch_type 3 0.5704 0.05143 1.4276 0.075 .
## trt:gender 1 0.1799 0.01622 1.3506 0.164
## gender:starch_type 3 0.6101 0.05501 1.5269 0.030 *
## Residual 66 8.7907 0.79263
## Total 75 11.0906 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There no significant difference in community composition and phylogeny between treatments but there is a difference between genders and at least 2 gender:starch_type pairs.
mouse_pcoa_uu <- ordinate(mouse, method="PCoA", distance=mouse_uu_dm)
plot_ordination(mouse, mouse_pcoa_uu, color = "trt") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal[3:4]) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/uu_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_uu, color = "starch_type") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal2) +
labs(col = "Starch Type") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/uu_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_uu, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/uu_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)there’s an interesting clustering, but we’d need more metadata to see what is causing this difference.
Weighted Unifrac
Merged
wu_dm <- phyloseq::distance(starch, method = "wunifrac")## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]adonis2(wu_dm ~ group*gender, data = data.frame(starchdf), permutations = 9999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = wu_dm ~ group * gender, data = data.frame(starchdf), permutations = 9999)
## Df SumOfSqs R2 F Pr(>F)
## group 3 4.0276 0.44175 30.339 0.0001 ***
## gender 1 0.5911 0.06483 13.357 0.0001 ***
## group:gender 3 0.3834 0.04205 2.888 0.0084 **
## Residual 93 4.1153 0.45137
## Total 100 9.1173 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis2(wu_dm ~ cohort*trt*gender, data = data.frame(starchdf), permutations = 9999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = wu_dm ~ cohort * trt * gender, data = data.frame(starchdf), permutations = 9999)
## Df SumOfSqs R2 F Pr(>F)
## cohort 1 3.9009 0.42786 88.1552 0.0001 ***
## trt 1 0.0558 0.00612 1.2619 0.2566
## gender 1 0.5791 0.06352 13.0872 0.0001 ***
## cohort:trt 1 0.0828 0.00908 1.8705 0.1291
## cohort:gender 1 0.1137 0.01247 2.5698 0.0636 .
## trt:gender 1 0.2141 0.02348 4.8382 0.0073 **
## cohort:trt:gender 1 0.0556 0.00610 1.2559 0.2478
## Residual 93 4.1153 0.45137
## Total 100 9.1173 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Significant difference in community structure and phylogeny between cohorts, genders and at least 2 trt:gender pairs.
pcoa_wu <- ordinate(starch, method="PCoA", distance=wu_dm)
plot_ordination(starch, pcoa_wu, color = "group") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal) +
labs(col = "Group") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/wu_group.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(starch, pcoa_wu, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/wu_gender.png", units = "px", width = w, height = h, create.dir = TRUE)Human
human_wu_dm <- distance(human, method = "wunifrac")## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]adonis2(human_wu_dm ~ trt*gender, data = data.frame(humandf), permutations = 9999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = human_wu_dm ~ trt * gender, data = data.frame(humandf), permutations = 9999)
## Df SumOfSqs R2 F Pr(>F)
## trt 1 0.02845 0.04848 1.2105 0.2718
## gender 1 0.03558 0.06063 1.5141 0.1968
## trt:gender 1 0.02929 0.04992 1.2465 0.2612
## Residual 21 0.49351 0.84097
## Total 24 0.58683 1.00000There is no significant difference in community structure and phylogeny between treatments, genders or the interaction between them.
human_pcoa_wu <- ordinate(human, method="PCoA", distance=human_wu_dm)
plot_ordination(human, human_pcoa_wu, color = "trt") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal[1:2]) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/wu_trt_human.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(human, human_pcoa_wu, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/wu_gender_human.png", units = "px", width = w, height = h, create.dir = TRUE)Mouse
mouse_wu_dm <- distance(mouse, method = "wunifrac")## Warning in matrix(tree$edge[order(tree$edge[, 1]), ][, 2], byrow = TRUE, : data
## length [3823] is not a sub-multiple or multiple of the number of rows [1912]adonis2(mouse_wu_dm ~ trt*gender*starch_type, data = data.frame(mousedf))## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = mouse_wu_dm ~ trt * gender * starch_type, data = data.frame(mousedf))
## Df SumOfSqs R2 F Pr(>F)
## trt 1 0.0982 0.02121 2.7985 0.027 *
## gender 1 0.6692 0.14455 19.0702 0.001 ***
## starch_type 3 0.6250 0.13501 5.9370 0.001 ***
## trt:gender 1 0.2444 0.05279 6.9650 0.001 ***
## gender:starch_type 3 0.6767 0.14616 6.4276 0.001 ***
## Residual 66 2.3161 0.50028
## Total 75 4.6296 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1pw <- pairwiseAdonis::pairwise.adonis2(mouse_wu_dm ~ starch_type, data = data.frame(mousedf))
print("significantly different pairs: ")## [1] "significantly different pairs: "for (i in range(2,(length(pw)))){
pair <- pw[[i]]
if (pair$`Pr(>F)`[1] < 0.05) {
print(names(pw)[i])
}
}## [1] "PTB_vs_BEP"
## [1] "CKP_vs_CTL"There is a significant difference in community structure and phylogeny between treatments, genders, at least 2 starch types, at least 2 trt:gender pairs, and at least 2 gender:starch_type pairs. The significantly different starch type pairs are: PTB-BEP and CKP-CTL.
mouse_pcoa_wu <- ordinate(mouse, method="PCoA", distance=mouse_wu_dm)
plot_ordination(mouse, mouse_pcoa_wu, color = "trt") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal[3:4]) +
labs(col = "Treatment") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/wu_trt_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_wu, color = "starch_type") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal2) +
labs(col = "Starch Type") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/wu_type_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)
plot_ordination(mouse, mouse_pcoa_wu, color = "gender") +
#stat_ellipse(alpha = 0.5) +
scale_color_manual(values = pal3) +
labs(col = "Gender") +
theme_bw() +
theme(
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
)
ggsave("../30_diversity-analysis/plots/beta-diversity/wu_gender_mouse.png", units = "px", width = w, height = h, create.dir = TRUE)